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Abstract 


This  paper  describes  a  method  of  obtaining  results  from  the 
simulation  of  a  countably  infinite  state  positive  recurrent  aperiodic 
Markov  chain  at  a  cost  considerably  below  the  cost  required  to  achieve 
the  same  accuracy  with  pure  random  sampling.  By  reorganizing  k  inde¬ 
pendent  epochs  or  tours  simulated  serially  into  k  replications  simu¬ 
lated  in  parallel,  one  can  induce  selected  joint  distributions  across 
replications  that  produce  the  cost-saving  benefits.  The  joint  distri¬ 
butions  follow  from  the  use  of  rotation  sampling,  a  special  case  of  the 
antithetic  variate  method.  The  chains  considered  are  of  the  band  type 
so  that  for  the  state  space  5  *  (0,1,2,...)  there  exists  an  integer  <5 
such  that  transition  from  a  state  i  can  move  no  further  than  to 
states  i  -  6  and  i  +  6  . 

The  paper  shows  that  an  estimator  of  interest  has  variance  bounded 
above  by  0(62(ln  k ) 4/k 2 )  when  using  rotation  sampling,  as  compared  to 
a  variance  0(l/k)  for  Independent  sampling.  Moreover,  the  mean  cost 
of  simulation  based  on  rotation  sampling  has  an  upper  bound  0 ( ( 61 n  k)2) 
as  compared  to  at  least  0(k)  for  independent  sampling. 

The  paper  also  describes  how  one  can  exploit  special  structure  in  a 
model  together  with  rotation  sampling  to  improve  the  bound  on  variance 
for  essentially  the  same  mean  cost. 

KEYWORDS:  Markov  chains;  rotation  sampling;  simulation;  variance  reduction. 
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Introduction 

A  recent  paper  (Fishman  1981)  describes  how  one  can  use  rotation 
sampling,  a  special  case  of  the  antithetic  variate  method,  to  induce 
substantial  variance  reduction  in  the  simulation  of  a  finite  state 
Markov  chain.  Since  many  discrete  event  simulations  have  an  underlying 
Markov  structure  or  one  close  to  being  Markov,  this  variance  reducing 
proposal  has  clear  appeal.  Moreover,  for  large  and  possibly  ill- 
conditioned  transition  matrices,  one  may  prefer  the  Monte  Carlo  or 
simulation  method  with  appropriate  variance  reducing  plans  to  numerical 
analysis  when  solving  for  steady-state  and  first  passage  time  distributions. 
In  fact,  it  may  be  the  only  feasible  method  for  some  problems.  The  present 
paper  extends  the  earlier  results  for  rotation  sampling  to  Markov  chains 
with  countably  infinite  state  space.  Here  the  willingness  to  rely  on  the 
simulation  method  for  solution  increases  because  the  convenience  of  an 
efficient  computer  code  for  matrix  inversion  is  no  longer  a  relevant  issue. 

The  earlier  work  derived  its  cost-saving  potential  from  viewing  the 
simulation  of  k  tours  in  series  of  a  finite  (n+l)  state  positive  re¬ 
current  aperiodic  Markov  chain  as  equivalent  to  the  simulation  of  k 
replications  of  the  Markov  chain  in  parallel.  Although  the  marginal 
distributions  that  arise  with  the  two  alternative  formulations  are 
necessarily  the  same  for  corresponding  variables,  the  parallel  formulation 
allows  one  to  induce  joint  distributions  across  replications  that  lead  to 
a  significant  cost  saving.  The  induced  joint  distributions  follow  from 


,ti 


A 


the  use  of  rotation  sampling,  as  described  in  detail  in  Fishman  and 
Huang  (1980).  The  cost  saving  arises  in  two  ways.  Firstly,  for  fixed 
n  ,  run  time  in  the  correlated  case  is  0(ln  k)  in  contrast  to  0 ( k ) 
for  the  serial  simulation.  Secondly,  for  fixed  n  the  variance  of  an 
estimator  of  interest  has  an  upper  bound  0((ln  k/k ) 2 )  for  the  corre¬ 
lated  case  compared  to  0(l/k)  for  the  serial  case. 

In  the  present  paper  we  replace  the  specification  n  <  »  with  the 
assumption  that  transition  from  a  state  can  go  to  no  more  than  26  +  1 
states  where  6  is  an  integer.  Then  it  is  shown  that  the  mean  cost  of  simula¬ 
tion  has  an  upper  bound  0((6  In  k)2)  and  the  variance  of  the  correspond¬ 
ing  estimator  has  an  upper  bound  0( 62 ( 1 n  k ) u/k 2 )  .  These  results  compare 
favorably  with  those  using  independent  sampling  to  simulate  the  behavior 
of  a  Markov  chain  where  mean  cost  is,  at  best,  0(k)  and  variance  is  proportion¬ 
al  to  1/k  . 

The  paper  also  shows  how  one  can  combine  rotation  sampling  with 
special  structure  in  the  chain  to  achieve  additional  variance  reduction 
without  any  essential  increase  in  cost.  The  relevance  of  this  result  be¬ 
comes  more  apparent  as  the  number  of  states  occupied  by  the  k  parallel 
replications  increases  at  a  given  step.  The  benefit  is  achieved  by  induc¬ 
ing  an  appropriate  joint  distribution  for  the  transition  paths  of  all 
replications  from  all  states  at  each  transition  while  preserving  the 
correct  marginal  distribution  for  each  replication.  By  contrast.  In  the 
earlier  use  of  rotation  sampling  we  preserved  the  marginal  distributions 
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but  merely  induced  an  appropriate  joint  distribution  for  the  transition 
paths  of  the  replications  in  a  given  state  on  a  given  transition,  leav¬ 
ing  the  sets  of  paths  for  different  exited  states  conditionally  indepen¬ 
dent  . 

Section  1  introduces  the  Markov  chain  notation.  It  also  formulates 
the  experiment  as  k  independent  tours  where  a  tour  begins  with  an  exit 
from  state  0  and  ends  upon  first  entry  into  state  a  .  A  reformulation 
in  terms  of  k  independent  parallel  tours  or  replications  is  presented 
and  then  extended  to  the  case  of  k  correlated  replications  using  rotation 
sampling.  Results  for  mean  number  of  transitions,  expected  cost  and  vari¬ 
ance  are  then  derived.  Section  2  derives  the  comparable  results  when  com¬ 
bining  rotation  sampling  with  special  structure.  Section  3  demonstrates 
how  the  results  apply  to  a  simulation  of  a  nearest  neighbor  Markov  chain. 

1 .  The  Infinite  State  Chain 

Consider  a  positive  recurrent  aperiodic  Markov  chain  with  state  space 
S  =  (0,1,2,...)  and  transition  probabilities  {p.^;  i , j  =  0,1,...} 
where  there  exists  a  positive  integer  6  such  that 

Pij  =0  for  |i-j|  >  6 

and 

i+6  ^ 

l  =  1  • 

j*max(0,i-6)  J 
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It  is  convenient  to  describe  an  alternative,  but  equivalent,  representa¬ 
tion  to  (1)  whose  value  is  apparent  when  actually  generating  sample 
paths  by  simulation  on  a  computer.  Let  Sj  denote  the  total  number  of 
states  that  have  positive  transition  probabilities  from  state  j  and 

let  m  ;  r  =  l,...,s.  denote  the  ordered  sequence 
J  J 

(m.  c  m .  , ;  r  =  1 , . . .  ,s  .  -  1 )  of  the  s.  states  to  which  entry  can 

J  *  J  »  ■  +  •  J  J 

occur  from  state  j  .  Then  one  has  the  representation 


d  .  >0 

•  im . 
j  jr 


r  =  1 . s  . 

J 


I 


s . 

J  P  • 
r=l  Jmjr 


=  1 


6  >max  (  Im^  -  j|  ,  |mjs  -  j|) 


j  »  0,1,... 


Let  A.,  denote  the  reward  received  when  a  jump  occurs  from  state 
^  J 

i  to  state  j  and  assume  for  the  moment  that  | A .  - 1  s  0(1)  .  Suppose 

*  J 

one  wishes  to  estimate  uQa  ,  the  mean  reward  received  on  a  path  begun 
with  a  departure  from  state  0  and  terminated  with  the  first  entry  into 
state  a  .  If  gives  the  number  of  transitions  from  i  to  j  on 

the  mth  of  k  independent  replications  and  S.  -  S  -  a  ,  then 

a 


«k-£i  i  i 

m=l  1 £Sa  jeS  0 

a 


(2) 


is  an  unbiased  estimator  of  Uga  and  var  «  1/k  .  Here  each 


A 


replication  begins  with  a  departure  from  state  0  and  ends  upon  entry 
into  state  a.  To  simulate  this  chain  on  a  computer  for  k  consecutive 
independent  replications  or  epochs,  one  need  only  set  p  n  =  1  and 
pa.  =0  for  j  =  1,2,...  . 

Let  denote  the  number  of  transitions  in  k  epochs.  Clearly 
E  Tk  =  0(k)  .  Let  denote  the  cost  of  simulating  k  replications 
in  series.  In  the  case  of  (1)  with  a  finite  state  space  of  n  +  1  , 
one  has  0(k)  <  E  S^  <  0(6k)  .  The  lower  bound  applies  if  one  can  store 
all  the  distributions  and  the  aliases  required  to  use  the  alias  method 
of  Walker  (1977).  Also  see  Kronmal  and  Peterson  (1979).  As  n  increases 
the  feasibility  of  this  approach  diminishes.  The  upper  bound  0(6k)  re¬ 
sults  from  application  of  the  inverse  transform  method  to  determine  the 
branch  taken  from  each  state.  Since  in  the  infinite  state  case  one  can 
conceive  of  using  the  alias  method  for  a  finite  set  of  commonly  entered 
states  and  the  inverse  transform  method  for  less  frequented  states,  taking 
E  Sk  >  0(k)  is  a  useful  bound  for  comparative  purposes. 

Consider  simulating  the  Markov  chain  with  replications  executed  in 
parallel.  Simulation  begins  with  an  exit  of  all  k  replications  from 
state  0  .  Then  one  sets  p_.  *  1  and  p  .  =  0  jeS.  ,  if  state  a  is 

3d  3  J  3 

to  be  absorbing  state.  Note  that  the  chain  is  absorbing  with  S  a 
transient  set  of  states.  The  simulation  ends  with  the  entry  of  the  last 
of  the  k  replications  into  state  a  .  Let  K,.^  denote  the  number  of 
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independent  replications  that  move  from  1  to  j  on  transition  £  and 
let  denote  the  number  of  replications  in  state  j  after  transition 
l .  Then 


l  l 


£  =  1  i.jeS 


AU 


is  an  unbiased  estimator  of 


■'Oa 


wi  th 


var  =  var 


Note  that 


l  l 


£ =1  i.jcS 


K. 


ij£ 


l  l 

m=l  IcS. 


(m) 


l  Nl 
j.S  ^ 


(3) 


Let 

Tk  =  min  (£:  K.^  =  k)  . 

From  Theorem  1  of  Fishman  (1981),  one  has  E  <  0(k)  .  Now  one 

easily  sees  that  at  transition  £  no  more  than  min(k,  <$£) 

states  are  occupied.  Therefore,  no  more  than  6t(t  +  l)/2  states 

are  occupied  in  t  transitions  and  the  expected  number  of  occupied 

states  has  upper  bound  *(E  f£  +  E  f^)  s  o(6k)  .  Let 

denote  the  cost  of  simulating  k  independent  replications  in  parallel 
and  note  that  no  more  than  26  +  1  states  can  be  entered  from  a  given 

state.  If  one  uses  a  sampling  program  as  in  Ahrens  and  Dieter  (1979) 

or  in  Fishman  (1979),  which  have  bounded  computation  times,  to  generate 


J 
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the  binomial  variates  K.  t  for  jeS,  r  =  and  =  1,2,..., 

Jmjr  J 

then 

E  Sk  .s  0(i?k)  . 

This  observation  is  peripheral  to  our  main  interest  and  we  merely  mention 
it  for  completeness . 

Once  the  replications  are  to  be  run  in  parallel  an  opportunity  ex¬ 
ists  to  induce  a  desirable  form  of  correlation  across  replications  while 
preserving  the  required  distributional  behavior  along  the  sanple  path  of 
each  replication.  Define 


V  *  -l., 


w  *  1 , . . .  ,s  .  ;  j  =  0,1 .. .  . 


(4) 


Let  U,U^,...,UK  be  i.i.d.  random  variables  from  u(0,l)  where 

j* 

KU  >  0  ls  givsn.  Then  for  parallel  replications  one  has 

r/  . 

o 

V’  “  V,  S.M.  V>(U“’ 


(5) 


where 


and 


"jO  -=  0 


U  <  X  <  V 


=  0 


f 


otherwise  . 
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Observe  that  on  each  transition  (£.)  a  Bernoulli  trial  determines 
the  path  that  each  replication  follows,  Independent  of  the  paths 


jm 


,lr 


v f  other  replications.  Moreover,  var(Kjm  ^(K^)  =  Kjj,  Pjm.  ^ 

Although  the  Bernoulli  property  must  be  maintained  to  assure  that  the 

sample  paths  for  each  replication  follows  the  correct  probability 

law,  the  independence  across  replications  is  but  a  consequence  of 

the  independence  o‘  l'  . U„  •  Hereafter  we  use  a  prime  superscript 

1  J* 

(')  to  denote  parallel  correlated  replications. 


)  . 


Suppose  that  for  r  3  1,...,K 

r-1 


/ 

J* 


Um  =  U  +  -K 


n 

+  2^1-1 
kTT 
j*. 


U  <  1 


r-1 


otherwise 


(6a) 


so  that 


4 

Sr1*1  '  S.r-r'hA’ 


(6b) 


We  refer  to  this  representation  as  rotation  sampling.  See  Fishman 
and  Huang  (1980).  Using  the  notation  x  =  x  (mod  1)  and  LxJ  =  x  -  x 
one  then  has  from  Theorem  5  in  Fishman  (1981): 
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This  expression  is  considerably  more  efficient  to  use  than  (6b)  is  in 
practice  for  its  cost  depends  only  on  the  number  of  occupied  states  and 
not  on  the  number  of  nonabsorbed  replications.  Hereafter,  we  assume 
that  (7b)  is  used. 

Let 


U 


-  I 


j  ,nuS 


A.  K#.  „ 
jm  jm£ 


/ 


(8a) 


(8b) 


/ 

Here  R^  denotes  the  sample  contribution  to  reward  on  transition  £  . 

/ 

Using  the  results  In  (7)  one  can  easily  show  that  var  Rk£  <;  0(6*£2), 
which  is  independent  of  k  .  In  the  case  of  a  finite  state  space 
Fishman  (1981),  shows  that  for  =  min(t:  Kgt  =  k)  and  cost  5^ 

ii  , 

E  Tk  s  ^  \  "  ®(ln  **)  ancl  var  ^k  5  0((lnk/k)2)  •  However,  one 
cannot  carry  over  these  results  directly  for  the  case  of  an  unbounded 
state  space.  Before  deriving  the  corresponding  results  for  this  case, 

t 

we  study  an  interesting  property  of  Tk  which  holds  for  both  finite 

and  infinite  state  models.  Because  of  the  restrictions  on  the  p..’s, 

•  J 

at  most  6  +  min(a,5)  transient  states  have  nonzero  entries  on  transition 

I 

Tk  -  1  .  In  fact,  absorption  can  occur  only  when  all  the  nonzero  entries 
are  in  states  a  -  min(a,6),...,a-l  and  a+l,...,a+6  In  the  way  specified 
i n  Lemma  1 . 
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Lemma  1 .  In  order  for  total  absorption  to  occur  on  transitions 
Jl+1 ,  4+2,...,  necessary  and  sufficient  conditions  at  the  end  of  transi¬ 
tion  l  are 

(a)  =  0  J j  -  a J  >6 

(b)  K'^<  1/(1  -  pja)  |j  -  a|  s  6  . 

Proof.  Condition  (a)  is  a  consequence  of  p.  =  0  for  | j - a |  >  6  . 
-  •  ja 

For  part  (b)  we  note  that  total  absorption  using  rotation  sampling 
implies 

Kja,i+1  =  LKj£  pjaJ  +  1  =  (j  -  a(  s  6  . 

For  a  specific  j  this  occurs  with  positive  probability  if  and  only  if 

i/n-Pj,)  • 

Lemma  1  provides  the  basis  for  an  initial  characterization  of 
absorption  time  in  Lemma  2. 

Lemma  2.  Let  k  >  1/(l-a)  where  a  =  sup  p2.  . 

Let 

t£  ■  min (t:  conditions  (a)  and  (b)  of  Lemma  1  obtain) 


and 
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i  * 

Yk  -  Tk  -  Tk 


Then 


(i) 

(ii) 


E  Tk  =  E  Tk  +  E  Yk  • 

is  independent  of  k. 


Proof.  Part  (i)  follows  by  inspection.  Part  (ii)  makes  use  of  the 

1  .  *  ★ 
fact  that  since  all  <  1/(1  -a)  k  for  i  =  +  1,  +  2,..., 

then  the  remaining  time  to  absorption  must  be  independent  of  k  . 

★ 

In  effect,  rotation  sampling  is  nonoperative  for  t  >  . 

Theorem  1  provides  a  more  comprehensive  characterization  of  absorption 
/ 

time  T, 


Theorem  1.  Let 


lo) 

p^j  =  probability  of  moving  from  i  to  j  in  i  steps. 


Mi  j  ft  =  Kfj2  "  Ki,£-1  Pij  ’ 


Then 


«>  4  ■  k  tof  +  Zjt 

where 


t-1 


Zit  =  l.  +  l  .  I  Hin,0  P, 


and 


*  W  ijt  Vl 

a  a 


4 •  k<>  - i  .  4l))  *  z 


(t-£) 


£=1 


Oa  '  at 


US, 


r 
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where 


'at 


=  l 


meS. 


mat 


t-1 

+  l  l 

£=1  i ,meS, 


Mim£ 


(1  -  l 


t-£ 


r=l 


p(r))  • 
Hma  ' 


(ii )  =  0(ln  k)  w.p.l . 

(iii)  E(T')r  =  0((ln  k)r)  r-1,2 . 


See  the  Appendix  for  proof. 

Theorem  2  shows  the  benefit  of  rotation  sampling  for  the  representa¬ 
tion  in  (1)  . 


Theorem  2.  For  the  simulation  of  a  Markov  chain  as  in  (  l  )  using 
rotation  sampling  as  in  (6a)  and  (7b)  : 

(i)  E  S'K  s  0(($ln  k)2)  . 

(ii)  var  Rk  s  0(62(ln  k)1*  /k2)  . 

E  S.  •  var  R, 

(m)  v(Rk«Rk)  =  e  $'  •  var  R£  *  0  ^2/(5  k ^  ^ 

t 

See  the  Appendix  for  the  proof.  The  expression  V(Rk,Rk)  gives  a 
measure  of  the  relative  efficiency  of  rotation  sampling  based  on  (7b) 
as  compared  to  independently  sampled  replications  for  simulating  the 
chain  (1)  . 


It  is  of  interest  to  compare  these  results  with  those  in  Fishman 
(1981)  for  an  n  +  1  state  chain.  Although  those  results  hold  for  the 
absorbing  state  a  =  0  ,  they  also  apply,  with  minor  adjustment,  more 
generally.  In  particular,  £  s  0(n2  In  k)  and 
var  R^  s  0((n  In  k)2/k2).  These  are  consistent  with  (1)  and  (11) 
of  Theorem  2  when  one  notes  that  the  mean  number  of  transients  states 
entered  has  an  upper  bound  min(0(n  In  k),  0(6(ln  k)2))  and  the  number 
of  states  to  which  a  transition  can  occur  has  an  upper  bound 
min(n+l ,  25+1 )  . 

/ 

Although  the  results  for  var  R,  hold  for  bounded  (A,.}  ,  the 

K  1 J 

/ 

results  for  the  variance  reduction  measure  V(Rk>  Rk)  and  the  absorption 
time  T^  apply  more  generally.  For  example,  the  relative  desirability 
of  the  proposed  sampling  plan  continues  to  grow  with  k  for  the  family  of 
reward  functions  (A..  = 

*  J 

reward  function  of  this  type  to  show  how  one  can  achieve  additional  variance 
reduction. 


r  p  ft 

cn  +  l  (c0  i  +  d.  j  )}  .  Section  2  uses  a 
u  Jt=l  *  * 


2 .  Exploiting  Special  Structure 

Although  Theorem  2  shows  that  for  large  k  rotation  sampling 
offers  a  clear  advantage  over  Independent  serial  replications, 
the  factors  6  and  (In  k)  in  the  bound  on  var  are  sufficiently 
broad  to  make  one  look  for  improved  convergence  for  moderate  k  . 

Consider  the  sampling  at  transition  t  .  Recall  that  rotation 

f 

sampling  applies  to  the  K  transitions  from  state  j  ,  but 

J  ^ 

transitions  from  different  states  are  Independent,  given  {K.s  ;  j  £s  }  . 

J  £  d 

We  now  describe  how  a  modification  of  this  independent  sampling  can 


i' 
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for  special  cases  of  (Ajm;  j,m  =  0,1,2,...}  lead  to  a  considerable 
reduction  in  variance.  Theorem  3  provides  the  basis  for  the  approach. 

Theorem  3.  Let 


0  <  p.|  <  1 

q_]“ 0 


qi  =  Vi  + 

i  ■  0,1, 

X’  ’  V,  ■  ^)(u> 

if  q,_,  i 

s 

if  q,_,  . 

Ys  =  l  Xi 

1=0  1 

s  -  1.2, 

where  U 

~  U(0,1 )  .  Then 

(1) 

YS  ■  LV  *  ,[0,qs)<U> 

5  ■  1.2,... 

(ii) 

Zs,t“W  tqtJ-  tqsJ 

where 

*  Js,t(U> 

‘  ICq5.qt)(U) 

if  s.  q 

if  qs  >  q. 

See  the  Appendix  for  the  proof.  The  significance  of  this  result 
becomes  apparent  shortly. 


(9) 
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Consider  a  Markov  chain  for  which  Sq  =  s^  =  . 


one  has 


Rk,£+1  =  ^j£Sa(Bj  Vj^+l  +  Aj  V} 


where  A.  =  A.  and  B.  =  A.  -  A.  for  j  =  0,1,. 
J  jm,,  J  jm.,  jm. 


n  ■  /•» 

j2 

B  1  :  0  .  Recall  that 


Jl  J  j2 


Sn^l+l  =  LKjaPjmjlJ  +  I[0,pj)(Uj)  * 


where  p.  =  K'a  pjm  and  Uq.U^,...  are  l.t.d.  from 

J 1 

that 


where 


Rk,£.+1  =  Cl  +  ij«Sa  8j  I[0.pj)(UJ) 


ct"  lusa  (BJ  Ll +  AJKJ*5  • 


=  max(j:  >  0)  . 


2  so  that 

(10) 

.  and 

11(0,1)  ,  so 

(11) 


Let 
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Then  from  Theorem  1 


j* 

J£ 


var(Rk,».llKji  J,sa>s7  T=06J  • 


(12) 


Now 


one  can  write  (VI)  equivalently  as 


RUi  ■  v  Vo (Bj ' Bj-'1  £~J  It0'p")(Un,) 


(13) 


where  one  should  note  that 


‘[o.pV'V  *  ‘tVl-OnA’ 


if  q  i  s  q 
m- 1  m 


*  1  -  *[q  ,q  ,)<Ua,>  ifVl>Rm 

LHm,Mm- r 


04) 


As  an  alternative  to  l.i.d.  Uq.U^,...,  set  Uq  e  U.j  = 
where  U  -  U(0,1)  .  From  Theorem  3  ,  one  has  (13)  as 


=  U 


iw  *  v  Vo(bj  ■  BJ-'!Civ 


irm  j  +  Jj-u;(u)]  •  "5> 


Now  observe  that 


V  £ 

,,r  <Rk,»l|Kjl  J'Sa>  s  I  V  IBj  ■  • 


(16) 
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which  is  to  be  compared  with  (12)  .  However,  a  more  appealing  case 

> 

exists.  If  for  j  *  0,1,...  either  Bj  &  BJ_1  or  s  Bj_1  ,  then 


one  has 


bj74  ■ 


As  an  example,  consider  the  reward 


Then 


Sr  ■  *  4  ‘2  V  *  c3 


BJ  *  AJ.j,  •  Ajmj2  •  C2(”jl  ‘  V 


SO 


that  (15)  has  for  either  m^  >  mj2  or  <  mj2  for  j  =  0,1,, 


var(Rk,t+l  lKj£  j£Sa)  5  tc2  (mjjl  "  tc2(Z<S  +  1)/2‘12  =  °(<5J) 


whereas  (11),  based  on  independent  UQ,U1 .  has 


c| 


var(Rk,«.+l lKjJt  jcV  *  T  Ij-0  (mjl  "  n,j2)2  *  +  l)s 
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The  extent  to  which  one  can  exploit  special  structure  for  more 

general  (Ajm  :  r=l . s^;  Sj=2  j=0,l,...}  and  for  the  case 

(s.  2  2}  remains  a  topic  for  future  research. 

J 

3.  An  Illustration 

This  section  describes  a  simulation  designed  to  show  how  the 
theoretical  results  of  Sections  1  and  2  fare  in  practice.  Consider  a 
single  server  queueing  system  with  independent  and  identically  distri¬ 
buted  exponential  interarrival  times  with  rate  X  ,  independent  and 
identically  distributed  exponential  service  times  with  rate  w  >  X  and 
infinite  capacity.  As  an  alternative  to  this  continuous  time  repre¬ 
sentation  one  can  view  this  system  as  a  nearest  neighbor  Markov  chain 

with  p.  .  ,  =  w/(X+u))  j  =  1,2,...  .  Here  5  =  1,  s.  =  2  for 
J  f  *  J 

j  =  1,2,...  and  the  absorbing  state  is  a  =  0  . 

The  objective  is  to  estimate  the  mean  number  of  customers  in  system, 
y  =  A/ ( to- X )  (e.g.  see  Gross  and  Harris  1974,  p.  67).  In  the  serial 

model  we  estimate  this  quantity  by 


A 


1 

,!  1  I  00  J(nW 

^  w  m»l  J»J+1 

1 

V 

A  m-l>l  '  J-J-' 

. 

0» 
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In  the  parallel  model  using  rotation  sampling  we  estimate  y 


by 


1 

IT 


= 


k  1  rk  r°°  ' 

A  +  A+w 


(18) 


Note  that  (17)  and  (18)  are  ratio  estimators  in  contrast  to  (2)  and 
(8b)  which  are  linear  estimators.  As  a  result,  (17)  and  (18)  are 
biased  estimators  of  y  .  Therefore,  our  evaluation  focuses  on  mean- 
square  error  rather  than  on  variance  alone.  The  motivation  for  consider 
ing  ratio  estimators  arises  from  the  observations  that  they  commonly 

arise  in  regenerative  simulation. 

For  convenience  and  without  loss  of  generality  we  set  w  =  1  . 

Table  1  gives  the  number  of  independent  macroreplications  performed 
for  each  value  of  k  and  experimental  layout  (p)  .  Table  2  presents 


Table  1 

Experimental  Layout 
(k  =  2m) 

Number  of  Macroreplications 


P  =  A/w 


1000 


100 


m  1,«.«,11 
m  —  !,«.«, 11 


m  *  12,. 
m  =  12,, 


.,14 

.,14 


f 


0.5 

0.9 


*  •  •  •  I 
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ratios  of  interest.  Here  c.  and  MSE^  denote  the  CPU  time  and  mean- 
square  error  for  model  i  where  i  =  1  denotes  k  independent  tours  in 
series,  i  =  2  denotes  parallel  replication  with  (7b)  and  i  =  3  denotes 
parallel  replication  using  (7b)  and  (14). 

For  p  =  0.5  Table  1  reveals  that  k  =  16  is. the  first  sample 

size  that  gives  benefits  for  both  models  2  and  3.  Also,  note  the 
strong  favorable  performance  for  k  2  256  .  For  p  =  0.9  the 
benefits  of  both  models  2  and  3  do  not  arise  until  k  =  1024  .  The  be¬ 
havior  thereafter  shows  substantive  variance  reduction.  In  presenting 
results  for  small  values  of  k  we  merely  sought  to  provide  a  comprehensive 
picture  of  how  the  models  behaved.  In  practice  one  normally  would  expect 
to  use  a  considerably  larger  k  so  that  the  risk  of  an  unfavorable 
variance  reduction  is  small.  Since  limited  time  was  spent  on  the  computer 
code  for  models  2  and  3,  one  also  suspects  that  careful  attention  to  pro¬ 
gramming  efficiency  would  move  the  indifference  points  to  smaller  values 
of  k  for  both  p  =  0.5  and  0.9. 


2 

4 

8 

16 

32 

64 

128 

266 

512 

1024 

2056 

4096 

8192 

16384 

32768 


(1) 

50 

58 

83 

1 

30 

1 

3? 

2 

23 

4 

53 

7 

16 

46 

37 

99 

71 

32 

230 

02 

451 

79 

696 

20 

2366 

17 

8526 

68 

.54 

.44 

.67 

1.00 

2.07 

3.69 

5.82 

14.71 

32.06 

64.38 

174.11 

583.63 

1476.74 

3783.70 

10280.72 


.22 

.23 

.17 

.22 

.29 

.33 

.42 

.54 

1.00 

1.95 

3.70 

4.98 

15.46 

17.42 


.22 

.20 

.19 

.28 

.25 

.33 

.42 

.67 

1.40 

2.07 

4.59 

10.38 

22.70 

42.82 

63.00 


65536 


22743.06 


37.21 


143.30 
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Appendix 

Proof  of  Theorem  1  (i)  comes  from  applying  rotation  sampling  to  k 

replications  in  state  j  =  0  on  transition  I  =  1  ,  then  applying  it  to 

each  for  states  j  =  0 . <$  on  transition  £  =  2,  etc.  One  then 

regroups  the  resulting  terms  to  give  {K*,  jcS)  . 

J  * 

For  (ii)  note  that 


usi  i‘S 


l  -  l  -  ■  0 


ieS  jeS 

a 


e  2jt  *  ° 


■k  "W  s  zjt  *  k<’  -  O 


JeS, 


and 


-w  - 1  .  pjf  > s  ht  5  k  l 


£=1 


£=1 


3U> 

30a 


Most  importantly,  observe  that  Z  ../k  does  not  depend  on  k  .  Also,  one 

J  ^ 


\  •  ■<"  2at  *  -k<’  '  It„  ■ 


has 
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Now  for  the  absorbing  state  there  exists  pc(0,l)  such  that 


L 


*  CXp1) 


so  that 


=  min[t:  In  0(0*")  =  In  (-Z^)  -  In  k] 


=  mi 


i n [ t :  t  In  0(p'^)  =  -In  (-Z^)  +  In  k] 


from  which  it  is  clear  that  =  0(ln  k)  w.p.l.  ( i i i )  is  a  direct 

/ 

consequence  of  (ii)  and  is  relevant  when  deriving  var  . 


Proof  of  Theorem  2.  Recall  that  no  more  than  6£  steps  are  visited  at 
transition  £  and  no  more  than  6t(t+l)/2  are  visited  in  t  transi¬ 
tions.  Also,  at  state  j  one  needs  to  determine  which  of  the  26  +  1 

f 

possible  states  are  entered  by  the  K.5  replications  at  transition 

J  ^ 

£  +  1  .  Using  (iii)  of  Theorem  1  gives 

E  s'  s  \  (26  +  1)6  E[T'(T'  +  1)  s  0( (6  In  k)2)  , 
which  proves  (i)  . 


For  (ii),  observe  that 


if'  -  k  n^-^n  +  U 

Kij£  "  k  P0i  pi j  ij* 


-27- 


where 


Wij£  =  Zi  ,2.-1  "  riji.  + 


r . 


ij* 


<k  Po!"”  4  Wij 


and  IK  is  the  uniform  deviate  used  to  determine  paths  from  state  i 
on  transition  i  .  Note  that  E  =  0  and  more  generally  that,  like 
Z.  ,,  W. .  is  not  a  function  of  the  magnitude  of  k  .  Recall  that 

I  i  Jc  I  •  J 

var  <,  0(62Z2 )  . 


Now,  one  can  write 


k£  =  ^  „  Aij  (k  P01  1)pij  +  WijJL} 


k£  .  *■  „  1j 

i«Sa  JcS 


and 


t  ,  #, 


var  Rk  .  E  '  v.r(Rk|T  )  ♦  E  •  Hf* 


k  k 


where 


V.-k  .  e(Vte> 


k  ieSa  jcS  £=1 


Since 


var(RkJljTk)  s  0( (6Tk>2) 


k|HT'|  <;  0(6(T.)2) 
'k  * 


and 
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one  has 


k2  var(R'|T')  s  0(52(t')4) 

k2H2'  £  0<62(t')4) 
Tk  K 


so  that 

k2  var  s  0(62(]n  k)4)  . 
Part  (iii)  follows  by  substitution. 


Proof  of  Theorem  3.  Clearly  Xi  ~  BeA(p^)  i  =  .  Let  r  =  L qs J 

Let  i i  denote  the  indices  at  which  q.  <  q.  .  for  t  =  l,...,r 

1 1  1 1“ 

Then  one  can  represent  Y$  as 


Ys  ■ 


C,  <ci  •  V.  5,t -,>  (un  *  C  KrV1”11 

1*1* 


Now  straightforward  evaluation  gives  (i)  so  that  Y$  -  lq$J~  Be^G^)  . 
In  particular,  var  Y$  =  qs(l  -  q&)  =  0(1)  . 

For  part  (ii)  we  have 


Y 


s 


Lqs2J  •  LqS]J  +  (U)  '  1  [0 ,<fs ^ )  (U) 
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